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We present an approximate metric for a binary black hole spacetime to construct initial data 
for numerical relativity. This metric is obtained by asymptotically matching a post-Newtonian 
metric for a binary system to a perturbed Schwarzschild metric for each hole. In the inner zone 
near each hole, the metric is given by the Schwarzschild solution plus a quadrupolar perturbation 
corresponding to an external tidal gravitational field. In the near zone , well outside each black hole 
but less than a reduced wavelength from the center of mass of the binary, the metric is given by a 
post-Newtonian expansion including the lowest-order deviations from flat spacetime. When the near 
zone overlaps each inner zone in a buffer zone, the post-Newtonian and perturbed Schwarzschild 
metrics can be asymptotically matched to each other. By demanding matching (over a 4-volume 
in the buffer zone) rather than patching (choosing a particular 2-surface in the buffer zone), we 
guarantee that the errors are small in all zones. The resulting piecewise metric is made formally 
C°° with smooth transition functions so as to obtain the finite extrinsic curvature of a 3-slice. In 
addition to the metric and extrinsic curvature, we present explicit results for the lapse and the shift, 
which can be used as initial data for numerical simulations. This initial data is not accurate all the 
way to the asymptotically flat ends inside each hole, and therefore must be used with evolution codes 
which employ black hole excision rather than puncture methods. This paper lays the foundations of 
a method that can be straightforwardly iterated to obtain initial data to higher perturbative order. 

PACS numbers: 04.25.Dm, 04.25.Nx, 04.30.Db, 95.30.Sf 


I. INTRODUCTION 

The simulation of binary black-hole systems is of fun¬ 
damental physical interest as the purely general relativis¬ 
tic two-body problem. It is also of astrophysical interest, 
since accurate simulations of the late inspiral and merger 
phases of such binaries will considerably help the effort 
to detect the gravitational-wave signals and extract infor¬ 
mation from them [J. Simulation reduces to the numer¬ 
ical solution of the Cauchy problem: take some initial 
data and evolve it. The evolution is difficult for many 
reasons, although in recent years there has been much 
progress. Still, any evolution is only as good as its initial 
data. 

A key issue of initial data is astrophysical realism. The 
goal is to compute data on a hypersurface that represents 
one moment in time of an astrophysical inspiral of two 
black holes. If such an inspiral is defined by initial condi¬ 
tions in the distant past for widely separated black holes, 
then the only way to obtain the exact data at a later time 
would be to perform the actual evolution using the full 
Einstein equations. This procedure, however, is compu¬ 
tationally expensive and thus impractical. On the other 
hand, several schemes have been developed to pose ini¬ 
tial data that approximates the astrophysical situation at 
a given time. These schemes are typically either based 
on post-Newtonian (PN) methods or on the numerical 
solution of the constraint equations of relativity. 

For example, the literature provides many types of ini¬ 
tial data for black holes in approximately circular or¬ 
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the constraints of the Einstein equations. To obtain such 
data, certain assumptions are made, such as conformal 
flatness and quasi-circularity. These assumptions are ex¬ 
pected to be good approximations within a certain error, 
although the astrophysical metric after a long inspiral is 
known to be not exactly conformally flat and the orbit 
not perfectly circular. 

In this paper we consider a post-Newtonian method 
combined with black hole perturbation theory to con¬ 
struct approximate inspiral initial data. For large to in¬ 
termediate separations of compact objects, an astrophys- 
ically relevant approximate spacetime can be obtained 
far from the black holes by analytical post-Newtonian 
and post-Minkowskian methods [l5|. One advantage of 
such methods is that they allow systematic improvements 
through higher order expansions (compared to numerical 
constraint solving schemes which typically include only 
the correct lowest order PN behavior). In their regime 
of validity, PN methods do result in appropriate devia¬ 
tions from conformal flatness and in non-circular inspiral 
orbits. The main disadvantage of PN methods is that, 
by construction, they are generally believed to fail in the 
final phase of the inspiral for fast moving objects, and 
also close to non-pointlike objects with horizons. On the 
other hand, black hole perturbation theory provides an 
accurate spacetime in a region of the manifold sufficiently 
close to the background black hole. The main disadvan¬ 
tage of this theory is that it fails sufficiently far from the 
background hole and, thus, cannot provide information 
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about the dynamics of the entire spacetime. 

In what follows, we take a concrete step toward com¬ 
bining PN theory with black hole perturbation theory 
using the mathematical machinery of asymptotic match¬ 
ing. The method maintains the potential for systematic 
improvement through higher order expansions, although 
we only work at low order here [la, [l7|. From the PN 
approach the method inherits its astrophysical justifica¬ 
tion, i.e. that for sufficient separation between the holes 
the method will yield metric components that are correct 
up to uncontrolled remainders of certain orders. The un¬ 
controlled remainders in the metric components have dif¬ 
ferent effects on different quantities, and it is of interest 
to see how other quantities such as the binding energy 
compare to those of other initial data sets in the litera¬ 
ture. This question is beyond the scope of this article, 
but should be addressed in the future. 

Concretely, we have to discuss how black holes are in¬ 
corporated in our approach. While formally PN meth¬ 
ods assume slow motion and weak internal gravity of 
the sources, it has been shown that the results hold 
as well for objects such as black holes with strong in¬ 
ternal gravity [l|[ as long as one is not too close to 
these objects. Near each black hole, a tidally perturbed 
Schwarzschild or Kerr spacetime provides another an¬ 
alytical approximation. Given that different approxi¬ 
mate metrics can be constructed from different scale ex¬ 
pansions, it is natural to try the method of matched 
asymptotic expansions [l9j | . If there is an overlap region 
(also known as a buffer zone) where both approximations 
(post-Newtonian and tidal perturbation) are valid, a dif- 
feomorphism can be constructed between charts used in 
different regions of the manifold by different approxima¬ 
tion schemes. Matching—demanding that both approx¬ 
imation schemes have the same asymptotic form in the 
overlap region—relates physical observables in the differ¬ 
ent regions, i.e. ensures that both expansions represent 
the same physical system. 

The first attempt to construct initial data in such a 
way was by Alvi [20}. By construction, there are discon¬ 
tinuities in the data, which were found to be too large for 
numerical experiments [2li |. Alvi’s fundamental problem 
was that, in the terminology of textbooks such as Ben¬ 
der’s [3, he did not match (construct expansions asymp¬ 
totic to each other everywhere in the overlap region) but 
rather patched (set approximate solutions to the Einstein 
equations equal to each other on specified 2-surfaces) so 
that large errors in the extrinsic curvature were possi¬ 
ble. Alvi’s Table I shows that his spatial metric near the 
black holes is discontinuous apart from the Minkowski 
terms (independent of G) in either region. Such discon¬ 
tinuities are problematic for numerical relativity, since 
part of the initial data is the extrinsic curvature which 
includes derivatives of the spacetime metric. Smooth¬ 
ing can be attempted, for example with transition func¬ 
tions as in Alvi’s next paper [22} . but it is not trivial to 
implement—especially with such large discontinuities— 
without adding unphysical content to the initial data. 


There is also the issue of making sure that the initial 
data slicing is treated consistently in the various expan¬ 
sions, which Alvi addressed to some extent but did not 
always make explicit his assumptions. Finally, there was 
a problem with the accuracy to which Alvi calculated 
metric components. Construction of the extrinsic cur¬ 
vature requires terms in the expansions that Alvi did 
not calculate because he assumed (incorrectly) that the 
counting of orders follows the standard pattern used in 
deriving post-Newtonian equations of motion. 

The main point of this paper is that we are able to 
correct the mathematical problems with Alvi’s approach, 
and that we provide initial data for actual numerical evo¬ 
lutions. We use true asymptotic matching to construct 
a piecewise metric for two black holes in circular orbit, 
including terms of order the gravitational constant G on 
the diagonal of the metric and 0(G ) 3 / 2 off the diago¬ 
nal. We then remove the piecewise nature of the ap¬ 
proximate metric by “merging” or “smoothing” the so¬ 
lutions in the buffer zones, thus generating a uniform 
approximate metric. We do this by constructing smooth 
transition functions so that the uniform approximate so¬ 
lution is in principle C°°, although in practice higher- 
order derivatives will be less accurate than lower-order 
ones. These transition functions are carefully constructed 
to avoid introducing errors in the smoothed global met¬ 
ric larger than those already contained in the approxi¬ 
mate solutions. This metric allows for the calculation of 
the lapse to 0(G), the shift to O(G) 3 / 2 , and the extrin¬ 
sic curvature to O(G) 3 / 2 . Although this data contains 
only the first order deviations from flat spacetime, our 
approach does include the tidal perturbations near the 
black holes. Strictly speaking, these tidal perturbations 
are valid only near the horizons—our approach cannot 
model the asymptotically flat ends inside the holes and 
therefore must be used in numerical evolutions with ex¬ 
cision rather than punctures. Our formalism can be ex¬ 
tended to higher order by including more precise post- 
Newtonian [23} and black hole perturbation theory re¬ 
sults [24| already in the literature. 

By construction, our initial data satisfies the con¬ 
straints only up to uncontrolled remainders of a certain 
order. Therefore, this data may still need to be projected 
to the constraint hypersurface via a conformal decom¬ 
position. One avenue worth exploring is that since the 
constraint violations are 0(G) 3 / 2 or smaller, it may be 
possible to find a constraint projection algorithm that 
changes the physical content of our initial data only at 
a comparably small order. In this manner, the formal¬ 
ism presented here can potentially be used to construct 
extremely accurate background data for constraint solv¬ 
ing. This hybrid combination of an accurate background 
4-metric plus constraint solving might potentially lead 
to very astrophysically realistic initial data, which then 
could be compared and tested against other sets already 
in the literature. 

This paper is organized as follows: Sec. [IT] describes 
the method of asymptotic matching as applicable to this 
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problem. Sec. EH discusses the near zone expansion of 
the metric and determines its asymptotic expansion in 
the overlap region. Sec. IIVI concentrates on the inner 
zone expansion of the metric and expands it asymptot¬ 
ically in the overlap region. Sec. [V] applies asymptotic 
matching to the metrics to obtain matching relations be¬ 
tween expansion coefficients and a map that relates the 
charts used in the different regions. Sec. |VI] constructs 
the global metric, discusses its properties, and builds 
transition functions to eliminate discontinuities between 
local approximations. Sec. IVIII computes the extrinsic 
curvature, lapse and shift. Sec. IVIIll concludes and points 
toward future research. 

Throughout this paper we use geometrized units (G = 
c = 1 ) and we have relied heavily on symbolic manipu¬ 
lation software, such as MAPLE and MATHEMATICA. 
We use the tilde as a relational symbol such that a ~ b 
means “a is asymptotic to 6 ” 0]. When we refer to our 
results as “global,” we mean that they cover the region of 
most interest to numerical relativity. Strictly speaking, 
our results do not cover the radiation zone (further from 
the binary than a reduced wavelength) or the asymptot¬ 
ically flat ends inside the holes. However, obtaining the 
radiation-zone solution and matching it to the near-zone 
solution is a solved problem 0 | and the asymptotic ends 
inside the holes can be removed with excision before nu¬ 
merically evolving the initial data. 


II. APPROXIMATION REGIONS AND 
PRECISION 

Let us now consider a binary black hole spacetime, 
with holes of mass mi and m 2 , total mass m = m 1 + 
m 2 and spatial coordinate separation b. The manifold 
(Fig. EJ) can be divided into 4 submanifolds, the bound¬ 
aries of which cannot be determined precisely due to the 
presence of uncontrolled remainders in black hole pertur¬ 
bation theory (BHPT) and post-Newtonian (PN) asymp¬ 
totic series. Nonetheless, an approximate subdivision is 
possible and we make one as follows: 

1. The inner zone of Black Hole 1, (submanifold 
Ci): f?i -C b, where Ra is the distance from the Ath 
black hole in isotropic coordinates. In this region, 
the metric is obtained via black hole perturbation 
theory as an expansion in e(i) = R\/b. 0,0. 

2. The inner zone of Black hole 2, (submanifold 

C 2 ): f ?2 6 , where the metric is obtained in the 

same manner as in region Ci but with labels 1 and 
2 swapped. 

3. The near zone, (submanifold C 3 ): va tua and 
r < A/ 2 - 7 T, where A is the wavelength of gravi¬ 
tational radiation, r is the distance from the bi¬ 
nary center of mass in harmonic coordinates, and 
ta — fnA is the separation in harmonic coordinates 
from the horizon of the Ath black hole. In this 



FIG. 1: Schematic diagram of the near zone (dark gray), in¬ 
ner zones (light gray) and buffer zones (checkered) projected 
onto the orbital plane. The black holes’ horizons are shown 
by solid black lines, while the excision boundaries are shown 
by dashed black lines. The near zone overlaps the inner zone 
of each black hole, and these overlap regions are the buffer 
zones (checkered patterns.) The boundaries of all zones are 
somewhat imprecise since they are based on power series ap¬ 
proximations, but the buffer zones are roughly spherical shells 
shown in this Figure as annuli. 

region, a post-Newtonian approximation is used 
for the metric with an expansion parameter e( 3 ) = 
via/ta 01 which is formally treated as the same 
order for both values of A. 

4. The far zone, (submanifold C 4 ): r > A/27T, where 
the metric is obtained from a post-Minkowski cal¬ 
culation [ 0 . 

These zones are shown schematically in Fig. |TJ projected 
onto the orbital plane. In these figures, the near zone is 
shown in dark gray, the inner zones in light gray and the 
buffer zones in a checkered pattern. The holes’ horizons 
are denoted by black solid lines, while the dashed black 
lines are the excision boundaries. 

The manifold is subdivided in such a way so that ap¬ 
proximate solutions to the Einstein equations can be 
found in each region. These approximate solutions will 
depend on certain coordinates and parameters local only 

(3) 

to that region. The near zone metric, gjj.J , is an expan¬ 
sion in £( 3 ) = m/r -C 1 , which depends on harmonic 
coordinates x M and parameters, such as the masses rriA 
of the holes and the angular velocity lo of the system. 
Similarly, the metric in inner zone 1 (or 2), /ij~ (or h ~), 
is an expansion in enj = r\/b <C 1 (or e( 2 ) = r-ilb <C 1 ), 
which depends on isotropic coordinates x& and certain 
parameters, such as the mass of the background hole M 
and the angular velocity Q of the tidal perturbation. The 
parameters and the coordinates used in different regions 
are not identical and are valid only inside their respective 
regions, although those regions overlap. 

A global metric can be obtained by relating the dif¬ 
ferent approximate solutions through asymptotic match¬ 
ing. The theory of matched asymptotic expansions was 
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first developed to perform multiple scale analysis on non¬ 
linear partial differential equations and to obtain global 
approximate solutions 0. In general relativity, this 
method was first applied by Burke and Thorne [H - 
Burke [ 27 j]; and D’Eath [HI, | 29 | in the 1970 s to derive 
corrections to the laws of motion due to coupling of the 
body’s motion to the geometry of the surrounding space- 
time. Based on these ideas, in this paper we will develop 
a version of the theory of matched asymptotics that is 
useful to obtain initial data for numerical relativity sim¬ 
ulations. 

Asymptotic matching consists of relating different ap¬ 
proximate solutions inside a common region of validity. 
This region is usually called the buffer zone by relativists, 
but is called the overlap region by others. For a binary 
there exists three such regions, which are 4 -volumes: Two 
buffer zones ( 0 13 and O23) are defined by the intersec¬ 
tion of the near zone and the inner zones of black hole 
1 and 2; the third one is defined by the intersection of 
the near zone and the radiation zone (O34). The former 
two, shown in Fig. |T]in a checkered pattern, are defined 
by the asymptotic condition uia ta b. The latter 
has been analyzed in Ref. [H and will not be discussed 
here. In this paper we perform asymptotic matching in 
the former two buffer zones 0 13 and 0 23 . In order for 
our tidal perturbations in the inner zones to be valid, the 
inner zones C\ and C 2 cannot overlap. 

Once a buffer zone has been found, asymptotic match¬ 
ing can be used to relate adjacent approximate solutions. 
The first step is to find the asymptotic expansions of the 
approximate solutions inside the buffer zones. These ap¬ 
proximate solutions depend on the expansion parameters, 
e (i)i e (2) and £(3), which are small only in their respec¬ 
tive regions of validity C±, C 2 , and C3. By definition, in 
each overlap region both expansion parameters are small, 
specifically em 1 and e( 3 ) <C 1 in O13, while e( 2 ) -C 1 
and 6(3) <C 1 in O2 3- Inside buffer zone 1 , for example, 
we can then asymptotically expand the near zone solution 

~( 3 ) 

in £(i) <C 1 to obtain gfiJ and the inner zone solution in 

6(3) <C 1 to obtain li~. These asymptotic expansions of 
approximate solutions are naturally bivariate since they 
depend on two independent expansion parameters. When 
working with these bivariate expansions, we use the sym¬ 
bol 0 {p , q) both to denote terms of order ( m/b) p (rA/b) q 
and uncontrolled remainders of order ( m/b) p or ( rA/b) q . 
Relating adjacent approximations then reduces to impos¬ 
ing the asymptotic matching condition 


g( 3 ) 


~ dx^ dx v 
dx^ dx v ' 


(1) 


This expression means not that the two approximate so¬ 
lutions are equated, which is correctly called “patching,” 
but rather that all coefficients of all controlled terms in 
the bivariate expansions are equated. 

After imposing the asymptotic matching condition, 
one obtains a coordinate and parameter transformation 
between the near zone and the inner zone 1 in Oi 3 (and 
similarly in 0 2 3.) These transformations allow for the 


construction of a global piecewise metric, which is guar¬ 
anteed to be asymptotically smooth in the buffer zone 
up to uncontrolled remainders in the matching scheme. 
Asymptotic smoothness here means that adjacent pieces 
of the piecewise global metric and all of their derivatives 
are asymptotic to each other inside the buffer zones. This 
asymptotic smoothness, however, does not rule out small 
discontinuities on the order of the uncontrolled remain¬ 
ders in the approximations, when we pass from one ap¬ 
proximation to the other. The global metric can be made 
formally C°° by smoothing over these discontinuities, 
which introduces a new error into the solution. Asymp¬ 
totic smoothness, however, guarantees that this error will 
be smaller than or equal to that already contained in the 
uncontrolled remainders of the approximations, provided 
the smoothing functions are sufficiently well-behaved. 

Finally, we enumerate the orders of approximation 
used in the near zone. The Einstein equations are guar¬ 
anteed to generate a well-posed initial value problem for 
globally hyperbolic spacetimes [l[ , where the initial data 
could consist of the extrinsic curvature Kij and the spa¬ 
tial 3 -metric 7 jj = gij. We can write the extrinsic curva¬ 
ture in the form 

Kij = — (27?(i/3j) — dtjij) , (2) 

where ( 3 i = goi is the shift vector and a is the lapse. 
Time derivatives are smaller than spatial derivatives by 
a characteristic velocity, which by the virial theorem is 
C^m/fr) 1 / 2 . Therefore to compute Kij consistently to 
a given order in m/6, the 4 -metric components goi are 
needed to 0(m/b j 1 / 2 beyond the highest order in g,j (and 
goo , which appears in a). In this paper we compute the 
first two nonzero contributions to the 3 -metric and ex¬ 
trinsic curvature, i.e. the leading order terms and the 
lowest-order corrections. This means that we need the 
4 -metric components goo and g^ to 0 (m/b ), but we need 
goi to 0 (to/&) 3 / 2 . Note that this does not correspond to 
any standard post-Newtonian order counting or nomen¬ 
clature, which is why we quote precisions and remainders 
precisely in terms of expansion parameters rather than in 
ambiguous terms such as “nth PN”. The standard post- 
Newtonian order counting corresponds to the calculation 
of the equations of motion for spinless bodies, but the 
counting must be altered when studying other problems, 
such as the bending of light or the equations of motion 
for bodies with spin [ 30 |. 


III. NEAR ZONE METRIC 

In this section, we present the post-Newtonian (PN) 
metric in the near zone C3 and expand it in the overlap 
region O13, the buffer zone where the inner zone C\ of 
BH 1 and the near zone overlap. When performing the 
matching in Sec. 0 we will obtain the corresponding 
expansion in the other overlap region O2 3 by a simple 
symmetry transformation. 
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We use harmonic corotating coordinates it, x, y, z) ro¬ 
tating around the center of mass such that 

t' = t , 

x' = x cos ut — y sin cut, 
y' = x sin u>t + y cos cot, 
z' = z, (3) 

where primed coordinates are nonrotating. The near¬ 
zone metric takes the form (3l| 


Implicit in this metric is the assumption that ta (for 
A = 1,2) is of order b, or in other words that the field 
point is not too close to one of the holes. (Recall that, 
in harmonic coordinates, the horizons are at ta = wa 
if we neglect tidal deformations.) We do not include 
the 0{m/r) 2 terms in goo in what is commonly called 
the first post-Newtonian (1PN) metric and we do include 
0{m/r ) 3,/2 terms in g 0 i for reasons discussed at the end 
of Sec. |TTJ In Eq. (|4|) it is sufficient to use the first post- 
Newtonian approximation 


(3) 

9oo 
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9 oi 
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( 4 ) 


where all remainders are at least 0{m/r) 2 and the su¬ 
perscript reminds us we are working on submanifold C 3 . 
Here 





( 6 ) 


to the angular velocity, where y = mim 2 /(mi + m 2 ) 
is the reduced mass of the binary. Our choice of sign 
corresponds to coordinates in which the orbital angular 
momentum is in the positive z direction, as shown in 
Fig. [2] 

We now concentrate on the overlap region (buffer zone) 
O 13 . Inside 0 13 we expand l/r 2 as a power series in ri 
as 


ri 


T2 


x\ + y 2 + z 2 = \J (x — m 2 &/m ) 2 + y 2 + z 2 , 

/x 2 +y 2 + z 2 = \/(x + mib/m) 2 + y 2 + z 2 (5) 




n —0 


( 7 ) 


are the usual harmonic radial coordinates centered on where the P n are Legendre polynomials. Substituting 
black holes 1 and 2 and b is the separation between holes. into Eq. 0 , we obtain 
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9 00 
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where all errors are of order 0(2,3) and where mi <C 
r\ <C b. The metric 0, denoted with a tilde, is the 
asymptotic expansion in the inner zone of BH1 of the 
PN metric, which is already an asymptotic expansion in 
the near zone. 


Observe that these expansions constitute a series 
within a series (bivariate series). In order to see this 


more clearly, we can rearrange the spatial metric to get 
2 mi 


~(3) 

9ij 


Sij | 

\ + 

Pi 

f Xl 


\n 


r 1 


1 _|_ 1712 n 
mi b . 


n 


'-T (9 > 


mi C n < b. 

Equation 0 is a generalized Frobenius series [32j, where 
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FIG. 2: Diagram of the near zone coordinate system used 
for the post-Newtonian metric in harmonic coordinates. The 
z axis is chosen parallel to the orbital angular momentum of 
the binary, so that the black holes orbit counter-clockwise here 
with angular velocity u>. The origin is chosen to be the center 
of mass, while n and denote the coordinate separations of 
the field point from hole 1 and 2 respectively. The coordinate 
separation between the holes is b. 

the expansion is about the regular singular points r± = 0 
and rq = oo. There are clearly two independent per¬ 
turbation parameters, namely 6 ( 3 ) = m\/r\ (the usual 
PN expansion parameter used in C 3 ) and e(i) = r\/b (a 
tidal perturbation parameter used in Coo). In the overlap 
region O 13 we can expand in both. 

IV. INNER ZONE METRIC 

In this section we discuss the metric in the inner zone 
Ci of BH1 and its asymptotic expansion in the overlap 
region O 13 . 

Physically, we expect the spacetime in the inner zone 
of BH1 to be Schwarzschild with mass Mi plus a tidal 
perturbation due to BH2. Thorne and Hartle [l 8 | argue 
that, in the local asymptotic rest frame (LARF) of BH1, 
the metric can be expanded in powers of M\ outside the 
horizon of BH1. The first term, independent of Mi, can 
be taken to represent the external universe and thus can 
be computed by placing a test particle in the spacetime of 


BH2 as done by Alvi [20| . This is the tidal perturbation 
due to BH2. Terms of higher order in M\ describe BH1 it¬ 
self (the Schwarzschild metric) and interactions between 
BH1 and BH2 (tidally-induced quadrupole, etc). At the 
level of approximation of this paper, we can neglect the 
interaction terms because they are 0 ( 2 , 1 ) or higher. 

Alvi identifies LARF coordinates, in terms of which 
the tidal perturbation is obtained, with isotropic coor¬ 
dinates (Fig. [3]). Observe that this coordinate system is 
centered on BH1 and is inertial. The asymptotic form of 
the inner-zone tidally perturbed metric valid in the buffer 
zone is given by Eq. (3.14) of Alvi [2C)|, who derives it by 
extending Thorne-Hartle type arguments. But to serve 
as initial data the perturbation is needed throughout all 



FIG. 3: Coordinate system used in inner zone 1. The isotropic 
coordinates are centered on black hole 1 and the other hole 
orbits with angular velocity Q. The matching parameters Mi 
and fl and the coordinates need not be equal to those 
used in the near zone. 


the inner zone including the strong-field region, not just 
in the buffer zone where the field is weak. Alvi derives 
a form of the perturbation valid everywhere in the inner 
zone and presents it in Eq. (3.23) of Ref. [2(3]. The result 
depends on parameters Mi and fl which will be related 
to the near-zone parameters mi and to when we perform 
the matching. 

We write Eq. (3.23) of Ref. [2C| as 
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( 10 ) 


where once more the superscript (1) is to remind us that 
this metric is valid in submanifold C 1 , while the IC super¬ 
script refers to isotropic coordinates. The uncontrolled 
remainders in all components of Eq. (ITT))) are 0{Ri/b) 3 . 

Inspection of Eq. m shows that the perturbed met¬ 
ric diverges as i?i —► 0 faster than (l/i?i) 4 , which pre¬ 
vents the use of puncture methods [lj|. Physically, this 
is because the tidal perturbation tacitly assumes a small 
spacelike separation from the event horizon. Of course 
this assumption is violated as i?i —> 00 , the asymptoti¬ 
cally flat spatial infinity we call the “outside” of hole 1; 
but it is also violated as Ri —* 0 since, in isotropic coordi¬ 
nates, that is the other asymptotically flat spatial infinity 
“inside” hole 1. Outside the hole we match to the near 
zone metric which is well behaved, but inside we have 
nothing to match to. Thus numerical evolutions using 
our initial data will need to excise the black holes. Ex¬ 
cision in practice still requires a slice that extends some¬ 
what inside each horizon, which raises the question of 
how far inside our data can be considered valid. Outside 
hole A the tidal perturbation is valid for Ra <C (rriA/m)b. 
The corresponding limit inside the hole is approximated 
by the transformation Ra —> (toa/2) 2 / Ra, which is an 
isometry for an unperturbed Schwarzschild hole, where 
we have used niA ~ Ma- Thus the tidal perturbation is 
good roughly for 


mA m ^ d „ , 

2 2 b m 


(ii) 


and the excision radius can be chosen anywhere between 
the lower limit of this expression and the horizon. 

Since matching will be simpler if performed between 
two coordinate systems that live in charts that are similar 
to each other, we choose to corotate first. We define inner 


isotropic corotating coordinates (ICC) 


X = X cos fiT + Y sin (AT, 

Y = —X sin fiT + y cos fiT, 

Z = Z, 

T = T. (12) 


Using these equations, we can obtain the inner metric in 
isotropic corotating coordinates, given by 


/l (1) * 

,L 00 

s H t + H s± n 2 (X 2 + Y 2 ) 


+2H st x2 (x 2 + F 2 - Z 2 ) , 

h w p 

rL ll n 

s H s 1 - H s2 - -j", 

h {1) p 

YX 

a -H„iYSl - H at -p-, 

h {1) p 
n 22 n 

Y 2 

» H s 1 - H s2 —, 

h {1) P 
rt 02 " 

a H s iXn + ^(X 2 ~Z 2 ), 

h {1) P 
ri 33 n 

z 2 

w TTsi — H s 2 , 

h {1) P 
ri 03 n 

. ,, YZ 
“ b 2 ’ 

h (1) P 
ri 12 " 

. „ 

h {1) P 
n 23 

. rr ZY 

h {1) P 
a l3 n 

. H zx 


( 13 ) 

















where we use the shorthand 


H st 

H s i 


H s2 




2 m 2 


1 - 


1 + 


1 + 


Mi 

2R[ 
Mi 


Mi 

2i?l 




2i?i 


-2 


1 + 


M 2 

If 


2i?i / 

1 - Mi/2i?i 
1 + Mi/2i?i 

RiJ 


4i?i 


+ 2 


7712 R 2 P 



4to 2 Mi / 
W?i P2 V 




(14) 


and the errors are still 0(Ri/b) 3 . In Eq. (fT51) we have 
dropped the superscript ICC in favor of (1), which refers 
to submanifold C\. 

By expanding Eq. m in powers of M\/R \, which is 
permissible in overlap region (H 13 , we obtain 


i (1) 

'00 


(i) 

01 

(i) 

02 

t (1) 

l 03 


+n 2 (x 2 + f 2 ) , 

2 to 2 [m „ 
b 3 \ b Z ’ 


1 + 


2Mi 2?ri2 




+ 


6 3_i?lP2 l^i 


X 


(15) 


where the errors are 0(2,3). Like Eq. (j 8 j) , /if is a bivari¬ 
ate expansion in both 6 (i) = l?i /6 (valid in the inner zone 
Ci) and 6 ( 3 ) = M 1 /R 1 (valid in the buffer zone O 13 ). In 
other words, it is the asymptotic expansion in the buffer 
zone to the asymptotic expansion in the inner zone. 


V. ASYMPTOTIC MATCHING 


In this section, we concentrate on finding a matching 
condition (^ 13 ) and a coordinate transformation (^> 13 ) 
that maps points in Ci labeled with isotropic corotat¬ 
ing coordinates (ICC) X M to points in C 3 labeled with 
harmonic corotating coordinates (HCC) x^. As already 
discussed, we concentrate on buffer zone O 13 , while the 
matching condition 023 and the coordinate transforma¬ 
tion </> 2 3 in (H 23 will be given later by a symmetry trans¬ 
formation. 


We assume that the coordinates are asymptotic to 
each other and that they can be expanded in an im¬ 
plicit bivariate series. That is, we assume that the map 
013 : X 1 ' —> x M has the form 


X 

Y 

Z 

T 


7/126 \ 
— 

m ) 




1 + 

( t . 

1/2 

) 7 i(z m ) + 

(?) 

1 + 

(¥, 

1/2 

Cl (*'*) + 

m 


' m 2 \ 

1/2 

7/1-2 \ 

1 + ( 

CT ) 

nix 11 ) + | 

— ) 


V b J 


(16) 


where y, 7 , 0 and r are functions of the harmonic corotat¬ 
ing coordinates which do not depend on m/b (or equiva¬ 
lently w), but are power series in r\/b. We continue these 
power series only to 0 (ri/ 6) 2 so that the errors here in 
013 are 0(2,3) as in Eqs. (0) and (fl5l) which are linked 


by 0 i 3 - The first term in the X equation above is chosen 
so that both coordinate systems have their origins at the 
center of mass of the binary. 

Like the coordinates, the matching parameters in the 
two coordinate systems must be identical to lowest order 
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in 7712 /b. (They must also be independent of coordinates.) 
Then rp\ 3 is given by 


Mi 

0 


777i 


U! 


( 7772 \ 1 / 2 , 7772 , / 777 2 \ 3 / 2 

[' + (—j ’» + X' B + (t) ’» 

( m 2 \ 1 / 2 , 777 2 , /777 2 \ 3/2 

1 + (—J K1 + ~r K2 + rrJ K3 


(17) 


(fl5l) . and (fl9l) at O^m/b) 1 / 2 are 


XI ~ 

-XXl,x, 

71 ~ 

-2/71.3/, 

tri iX 1 

~ xXi,t, 

tr^ z 

~ <1.7, 

>l,a: ^ 

-XXl , 2 , 


Cl-<1 ,z 

t\ ~ 

2 /71.t - tT h y 


/ 777-2 \ J / 2 

V m ) ’ 


®Xi,y ~ -2/7i,x, 

2/7i,2-<1,3/. (20) 


where the errors are 0 (rn/b) 2 . 

Before moving on with the calculation, let us discuss 
the physical meaning of the assumptions we have just 
made. Equation m implies that inner and near zone 
metrics are identical in the buffer zone up to a change in 
coordinates. For a single black hole, in the buffer zone 
(which is outside the horizon), the only difference be¬ 
tween isotropic and harmonic coordinates is the radial 
transformation 


/ M 2 \ 

r = R ( 1 + J^)’ (I8) 

where r is in harmonic coordinates and R is in isotropic 
coordinates. This has the asymptotic form posited in 
Eq. (fl6l) . Thus the assumption of Eq. (flGl) is only needed 
in the buffer zone for the matching in this section. How¬ 
ever, we will want to write our final results in a global 
coordinate system which goes inside the horizons (though 
not all the way to the asymptotically flat ends). For this 
purpose we assume that the form of Eq. m holds for 
all values of n > 0. This has the effect of defining a 
new coordinate system which is asymptotic to harmonic 
coordinates in the near zone and to isotropic coordinates 
in the inner zone. 

Now let us return to the asymptotic matching. Using 
Eq. (IT7I) we can transform Eq. (fl5l) to harmonic coro¬ 
tating coordinates and impose the matching condition of 
Eq. 0), 


where commas stand for partial differentiation. The so¬ 
lution in terms of integration constants Ci is 

tri(x,y,z,t) = C 4 x + C 5 y - \Jm 2 /my + C 8 z + C 9 , 
x\i (x,y,z,t) = -Ciy + C 2 z + C 4 t + C 3 , 

yii{x,y,z,t) = -Cix + C 7 z + C 5 t + C 6 , 

zCi(x,y,z,t) = -C- 2 x - C 7 y + C$t + C w . ( 21 ) 


F or simp licity, we choose all Ci = 0 except C 5 = 
\ J 7772/777. Thus the coordinate systems are identical at 
0 (tt 7 /&) 1 / 2 . The coordinate transformation then becomes 


X 

Y 

Z 

f 


m 2 b 


~LV b ) 

/777 2 \ 3 / 2 . 

(—) *,(*") 




1 +(^) 7 ^) + (^) 3 <x) 


1 7712 m 2 

b V m 


*+( 


^) 6 M + (^) 3/2 «*') 

1 + (x ) T+ (x )' T3 ^ 


( 22 ) 


where the errors are still 0(2,3). 

Applying asymptotic matching to 0(m/b), we obtain 


g( 3 ) 

fj IS 


(<) 




dX a dX? 

dxv dx v 


(19) 


Equation (1191) provides 10 independent asymptotic rela¬ 
tions per order, all of which must be satisfied simulta¬ 
neously. Each asymptotic relation results in a first-order 
partial differential equation for the coordinate transfor¬ 
mation, leading to 10 integration constants per order. As 
we shall see, these constants correspond to boosts, rota¬ 
tions, and translations of the origin. 

Equation m must be solved iteratively in orders of 
(■777/6 ) 1//2 . Evaluating the nonzero components (the di¬ 
agonals) of Eq. ([TUI) at zeroth order in m/b, i.e. com¬ 
paring Eqs. © and (El), provides no information, since 
it only asserts that at lowest order both metrics repre¬ 
sent Minkowski spacetime. This is true for any matching 
formulation involving metrics of objects that would have 
asymptotically flat spacetimes in isolation. 

The asymptotic relations given by evaluating Eqs. ©, 


- (t 2 + tT 2 , t ) 
X2 + XX2,x 
72 + 2/72.3/ 

C 2 + < 2,2 

XX 2 ,t - tT 2 ,x 


1 (x — vn 2 b /rri) _ 1 Xi 
b ~~ 6"’ 

1 (x — TO 26 / 777 ) 1 XI 

b = ~b’ 

1 (x — TO 26 / 777 ) 1 Xi 

b = ~b’ 

(x — 77726 /m) Xi 

b _ fe"’ 


- tT2,y '• 

\7772 

< 2 ,t '• 

- tT 2 , z , 

X\2,v - 

- -2/72,x, 

XX2,z '• 

^ -<2.*, 

72,2 B 

" -<2.y 


\ 1/2 

777 \ y 

- T Kl > 

7772 / O 

/ 777 \ 1 ^ 2 x — m 2 b /777 


-« 1 , 


( 23 ) 
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Once more we have a system of 10 coupled partial dif- = 0, so that the fi = ui to this order. The solution to 

ferential equations that now depends on «q, which is a this system is then given by 
parameter that relates fi and lo. For simplicity, we choose 

I 


t 2 = 


1 _ x m 2 
b m. 
x m 2 


+ Db l + ° 4 1 + Ai 7 + T’ 


l 


* 2 - + “ + ( 2i2 + y 2 + z 2 ) +d i- + d a - + d 2 - + 


2 xb 


z 

‘ X 


D, 


, x m 2 r. x „ z t D 6 

72 = 1 —— H- D \—h Dj —b D 5 —I-, 

b m y y y y 


C2 = 1- 


x m 2 


m 


z z z z 


(24) 


where the Di are 10 more integration constants. For simplicity we set them to zero, and the coordinate transformation 
becomes 


X 

Y 

Z 

f 


m 2 b\ 

-I ■+ 

m J 
(?) 
(?) 


‘M?) 


(?) 


1 - 


(x — 2 m 2 b/m) 
2b 


; +( ? r ^,} +i i (2 , + ^,), 


2/j 1 + 

z{ 1 + 


1 (x — m 2 b)/m 
b 

+ (?) s, w>] 

{x — m 2 b)/m 
b 

??) 3, ?<->j 

(x — m 2 b/m ) 
b 

??r^>j 



(25) 


with errors of 0(2,3). 

Now that matching has been completed to 0(m/b)°, 0{m/b ) 1 / 2 and 0(m/b), we can proceed with matching at 
0(m/&) 3 / 2 . However, keeping in mind our discussion in Sec. ns of the orders needed, we will only use the spatial- 
temporal part of the asymptotic relations matrix m, 


XX3,t - tr 3iX 

tT 3 ,y ~ 2/73,t 


- (—) 
\m 2 ) 


1/2 


4x 


m 2 


7 1 fc 2 - - -b 3- 


m 
m 2 

(? 


r X 

1 

2 ' 


1 

m 2 
m 


4mi - 6 m 2 / m 2 y 


2m i 


m 


V m 7 


+ 


+ 


(7 


m 2 

m 


m 

2m i 
m 


(!) 


*(i 

b 


4mi 2 m 2 

-b- 2 

m m 

4 ^b 


m 2 r i 


trs, z - zC 3,t 


777-2 

777- 


1 — AC 2 H - 2 


/m2 V 
V m ) 


m \ 1 ^ 2 zy 
m 2 y 6 2 


3m 2 — 4 fi 


9 p, + 8mi 

2m 


(26) 


For simplicity, we choose k 2 = 1 + 3 m 2 /m. This com¬ 
pletes the derivation of the matching parameters, since rji 
and «3 did not appear in the differential equations at all, 
and hence, we can neglect them to this order. Note that 
this choice of parameter matching is different from Alvi’s 
choice, and thus our coordinate transformation is also 
different. Up to 0(m/b) 2 , the corresponding parameter 


matching condition 1/73 is 


fi 


; U) 


1 + 


m 2 /l + 3 m2 

0 V 777- 


777- / J ’ 


Mi 


777i- 


(27) 


This choice of ipi 3 simplifies Eq. (1261) , which now becomes 
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tv.3,x - XX3,t 


tr 3 ,z - zC, 3>t 
tT 3 ,y - 2/73,t 



(28) 


As before, we choose the integration constants for simplicity (and to keep the slicings close to each other), resulting 
in the following solution: 


T 3 ( x , y , z , t ) 

X3{x,y,z,t) 

C 3 {x,y,z,t) 


73 (x,y,z,t) 



( 

r\ 2 (7 

4m 1 

2 m 9 \ 

+ - 

r - - 

-+- 

V 

bJ \2 

m 

m J 


m2 

m 


2m i 
m 


3(T) 2 + — 

3 \bJ m 2 ri 



9/i + 8m i 
2m 



(29) 


To summarize, we have found a coordinate transforma¬ 
tion <f> 13 and a set of parameter relations ^13 that pro¬ 
duce asymptotic matching to 0(3/2,3) in the 00 and ij 
components of the 4-metric and to 0(2, 3) in the Oi com¬ 
ponents. 

Note however that the 73 piece of the coordinate trans¬ 
formation becomes singular at r 1 = 0. Also recall that 
the point ri = 0 is not identical to the point Ri = 0, 
where the inner zone metric perturbation diverges. Hence 
if we excise the inner zone metric close to R\ = 0, the 
point r i = 0 might be outside the excised region, in which 
case our coordinate transformation would introduce a co¬ 


ordinate singularity outside the excised region. To get rid 
of this singularity we will replace ri by 

ri = ^ r\ + 6 m 2 . (30) 

This change amounts to adding a higher order term to 
the coordinate transformation, which has no effect in the 
buffer zone at the current order of approximation, but it 
has the advantage that the resulting coordinate transfor¬ 
mation is now regular at ri =0. With this replacement 
the coordinate transformation is given by 



12 


X 

Y 


Z 

T 

n 

M 1 




9y + 8mi | 3 
+ 2 


2 m 
z{ 1 + 


(x) 


^-(=?) 


1 - 


1 - 



1 + 


^ (l + 3—) 
b \ m J J 


TOi- 



(31) 


The coordinate transformation for matching in the other 
overlap region 023 is obtained by the following symmetry 
transformation: 

1 <-> 2, x —► —x, y —> —y, z z. (32) 

In Eq. PIT) , t should be considered small just as x, 
y : and z are. Recall that fundamentally the overlap re¬ 
gions are 4-volumes, although when we choose a time 
slicing we have to deal with their projections on a spa¬ 
tial hypersurface {itia -C r a b). Just as the over¬ 

lap regions span a limited range of xa, so they span a 
limited range of t. The post-Newtonian metric and the 
perturbed Schwarzschild metric are formally stationary, 
but the true physical system includes gravitational waves 
(not modeled here) which for example change the orbital 
frequency on a radiation reaction timescale. While this 
timescale is longer than an orbital period, which must 
be of order 6, rotation and boosts mix space and time 
terms and to be consistent with xa b we must keep 


t <C b. Therefore we choose the t = 0 = T slice when 
discussing the approximate metric in the next section, 
which restricts our overlap region to the intersection of 
this 3-surface with the overlap 4-volume. 


VI. AN APPROXIMATE METRIC FOR 
BINARY BLACK HOLES 

In this section we transform the inner zone metric in 
isotropic corotating coordinates to harmonic corotating 
coordinates via Eq. m- The metric in the inner zone of 
black hole 1 is given by 

g$ = hfh\j? v , (33) 

where in the buffer zone the goo and components have 
errors of 0(2,3) and goi has errors of 0(5/2, 3). In the 
above equation the Jacobian = d v X& can be ex¬ 
panded as 
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(34) 


Furthermore, Jijli, refers to the inner metric presented in 
Eq. (1T51) , where we substitute the coordinate transforma¬ 
tion given by Eq. m and 


R 1 = (X 2 +Y 2 + Z 2 ) 1/2 . (35) 


The metric in the inner zone of black hole 2 (C 2 ) is given 
by the symmetry transformation (|32D applied to Eq. 

We now have all the ingredients to construct an ap¬ 
proximate piecewise metric, for example 


piece 
y/iv 


o (1) 

7*1 

< rf 

(2) 

yX 

7*2 

< rf 

(3) 

yX 

7*1 

> rf 


(36) 

r < \/2it 


for some transition radii 7/4 which are chosen to be in¬ 
side Oa 3 for A = {1,2}. In Eq. (1551) . gjiJ is given in 
Eq. dH), g$ is given in Eq. (1551) and gj±J is the symmetry 


transformed version of Eq. (1331) . The inner zone pieces 
of this metric are accurate to 0(i?i/6) 2 , while the post- 
Newtonian near zone pieces are accurate to 0(rn/r) 3 ^ 2 . 
In the overlap region, the inner and near zone 0* compo¬ 
nents asymptotically match up to 0(3/2,2) as required 
for the extrinsic curvature, while all other components 
asymptotically match only up to 0(1,2). 

Note that when we applied the coordinate transfor¬ 
mation of Eq. m to the inner zone metric we kept 
terms up to 0(777/6) 5 , which at first glance seems too 
high. These terms are needed because the inner zone 
metric represents a tidally perturbed black hole with er¬ 
rors in the physics of O[(777/6) (i?i/6) 3 ]. Close to BH1 
(i?i ~ 777), the error in the physics is only of order 
0(m/b) A and hence we should keep terms to at least or¬ 
der 0{m/b) 3 . However, recall that the perturbed black 
hole metric satisfies the Einstein equations up to errors 
of only 0[(777/6) 5 / 2 (i?i/6) 2 ] even though its astrophysical 
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resemblance to a binary black hole has errors already at 
0[(m/b)(Ri/b) 3 ]. If we want to obtain a metric which 
is close to the constraint hypersurface, we should keep 
terms larger than 0[(m/b)(Ri/b) 3 ], but not larger than 
0[(m/b) 5 / 2 (Ri/b) 2 }. In particular, close to the black hole 
we have constraint violations of 0(m/b) 9 ^ 2 . For exam¬ 
ple, if we had dropped terms of order 0(m/&) 4 in the 
inner zone metric we would have introduced additional 
constraint violations at this order. 


A. Global character of the asymptotic metric 

In this subsection we plot Eq. m to describe some 
features of the approximate piecewise metric. We choose 
a system of equal-mass black holes mi = m2 = m/2 sep¬ 
arated by b = 10m, so that both holes are located on the 
x axis with BH1 at x « 5 m and BH2 at x w —5m. Fig- 
ures|l]and[5]show the xx and 00 components of the metric 
along the x axis for this system. (Other components of 
the piecewise 4-metric exhibit similar behavior.) In all 
plots we use a dashed line to denote the near zone metric 
{(jj'i }) and dotted lines to denote the inner zone metrics 
[gjA ') We choose the separation b = 10 m because it is 
near the minimum for which our formalism makes sense, 
and we plot on the x-axis because it is where the worst 
behavior occurs. The idea is to (i) show some practi¬ 
cal features of matching which have not been presented 
in the literature and (ii) demonstrate the limits of the 
method, particularly regarding the minimum separation. 

In these plots, we also include error bars that esti¬ 
mate the uncontrolled remainders in the approximations. 
These remainders can be approximated by 

2[(mi/n) 2 + (m 2 /r 2 ) 2 ] (37) 

in the near zone, 

2(m 2 /b)(R 1 /bf (38) 

in inner zone 1, and the same with 1 <-> 2 for inner zone 
2. The error in the near zone metric was estimated by 
looking at the next order (2PN) term in the metric com¬ 
ponents [3l| . [That term is much more complicated than 
Eq. (1571) , but is numerically about equal to it in the region 
plotted.] The error in the inner zone metric comes about 
because the next tidal correction for a single black hole 
of mass m in perturbation theory will be roughly propor¬ 
tional to ~Riojo,k- The error bars of the approximations 
are position-dependent, and thus provide a useful sign 
of where each approximation is breaking down. The PN 
error bars are larger near the holes than far away, and 
the BHPT error bars exhibit the opposite behavior. The 
error bars also provide an indicator of where both ap¬ 
proximations are comparably good: Neglecting angular 
factors (as is typical in the literature), the error bars for 
the near zone and inner zone A are comparable at 

r T A « (b 4 m 2 A /m) 1/5 , (39) 


which takes a value of about 4.8m for the system plotted 
here. This radius is a good candidate for the “transition 
radius” of Eq. ([36H . but note that, in principle, there is 
an infinite number of possible candidates, as long as they 
are in the buffer zone. Furthermore, note that this is 
not a “matching radius,” since there is no such thing. 
Matching asymptotic expansions, as opposed to patch¬ 
ing them as done by Alvi ppj . does not happen at any 
particular place in the buffer zone. Rather, it makes two 
expansions comparable throughout the buffer zone up to 
the uncontrolled remainders. 

In Figs. 31 and [5] we plot the xx and 00 metric com¬ 
ponents along the x-axis for the PN approximation as 
well as for the two perturbed black hole approximations. 
In Fig. □ the buffer zones around each black hole were 
sketched as spherical shells around the holes, formally 
defined by mi C ri C i and m 2 « r 2 < t. It is im¬ 
portant to recall that this definition, which is ubiquitous 
in the literature, is imprecise because of the <C symbols 
and one cannot simply substitute < symbols. (For one 
thing, there is angular dependence of the uncontrolled 
remainders.) Inserting the parameters of our system into 
these definitions, in Fig. [4] the intersection of buffer zone 
1 with the x-axis is given by 5.5 <C x/m -C 15 (to the 
right of BH1) and — 5 <C x/m <C 4.5 (to the left, be¬ 
tween the holes). In the part of the buffer zone to the 
right, away from BH2, we see a clean example of the be¬ 
havior expected of matched asymptotic expansions: The 
near-zone curve and the BHl inner-zone curve do not in¬ 
tersect, but the difference between them is comparable to 
the estimated error bars everywhere within this part of 
the buffer zone, even if we replace the -C operator in the 
definition of the buffer zone by the precise < operator. 

Between the holes, to the left in Fig. 31 the interpreta¬ 
tion of the curves is not so simple. We cannot replace -C 
by < in the definition of buffer zone 1 because to the left 
of x/m ~ 0 (ri ss b/2) the near-zone metric component 
resembles that of BH2 rather than BHl. This is because 
mi <C r± <C b is a rough criterion obtained by ignoring 
(among other things) the angular dependence of the ex¬ 
pansion coefficients in Eq. © , inspection of which shows 
about a factor of two variation as the angle is changed. 
This angular dependence means that if one tries to rede¬ 
fine the buffer zone heuristically as “the region where the 
error bars on two curves overlap and are not too large,” 
it is significantly aspherical and can be squeezed out en¬ 
tirely. Even at the origin (where they are smallest), the 
error bars from the PN approximation in the near zone 
are visibly larger than they are for most of the right- 
hand part of buffer zone 1. These error bars, however, 
are what is expected: At the origin in Fig. 3] the 0(m/r) 
term which is kept in the PN metric has a value of 0.4, for 
a total metric component of gxx = 1.4. The uncontrolled 
0(m/r) 2 term we use for the error bar is 0.04, which is 
precisely 10% of the m/b correction and about twice the 
value at x/m = ±10, a comparable distance on the other 
side of each hole. The fact that the error bars are worst 
in between the holes does not depend on m/b or the mass 
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FIG. 4: This figure shows the xx component of the near zone 
metric (PN), denoted by a dashed line, and the inner zone 
metrics (BHPT), denoted by dotted lines, along the (har¬ 
monic) x axis for a perturbative parameter m/b = 1/10, with 
the black holes located at x/m = ±5. The buffer zones cannot 
be precisely defined, but most of the region plotted is within 
one or the other (see text). Matching does not guarantee that 
two curves which are asymptotically matched intersect any¬ 
where in the buffer zone, but rather that they are comparable 
at the level of the uncontrolled remainders. The error bars 
estimate these remainders as described in the text. 


ratio, but rather is a reflection of the physical assump¬ 
tions on which matching is based. The near-zone metric 
is matched in buffer zone 1 to the metric of inner zone 1, 
and in buffer zone 2 it is matched to the metric of inner 
zone 2, but inner zone 1 is never matched to inner zone 
2. It is the intervening near zone that ensures that each 
black hole’s tidal perturbation is the appropriate one (up 
to uncontrolled remainders) for the other black hole, be¬ 
cause each tidal perturbation is derived for a black hole 
without a nearby body. (See, for example, the discussion 
in Sec. II.B of Thorne and Hartle fl8|.) 

Since the inner-zone and near-zone metric components 
diverge as va —> 0, Fig. 0] might seem to imply that the 
solutions approach each other near the horizons. This 
misconception can be rectified by scaling the solutions to 
the Brill-Lindquist factor ip 4 , where 

”= 1+ £: + £- < 4 °) 

This removes most of the divergent behavior of the so¬ 
lutions, as shown in Fig. [Gl In this figure, we only plot 
the region near BH 1 to show the difference in divergence 
better, but the region near BH 2 is very similar. To the 
right, away from the other black hole, we see that the 
near-zone and inner-zone solutions are indeed quite sim¬ 
ilar near BH1 and that there is a wide region where both 
sets of error bars are comparable and overlap. The tran¬ 
sition radius rf ss 4.8m. ( x/m ~ 10) discussed above is 
seen to be a good approximation of where the error bars 
are equal. To the left, between the holes, there is only 


FIG. 5: This figure is similar to Fig. [4] but it shows the 
00 component of the metric. Observe that in this compo¬ 
nent the differences between the different approximations are 
more pronounced, although the general features of asymptotic 
matching are still discernible. 
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FIG. 6: In this figure we plot the xx component of the near 
zone (PN - dashed line) and inner zone 1 (BHPT -dotted 
line) metrics divided by the Brill-Lindquist factor ip 4 . The 
behavior of the solutions is clearly different as we approach 
the event horizon from the left (the direction of the other 
hole). 
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a small region (about \x\/m < 1) where the error bars 
overlap, and they are never equal. 

The smallness of the region where the error bars over¬ 
lap is an indication that b = 10m is approaching the min¬ 
imum separation for which our approximation method 
makes sense. The disappearance of such a region could 
be used as a criterion for the failure of matching, al¬ 
though this is not a standard test and several different 
approximate criteria could be used (and this region is not 
the formal definition of a buffer zone anyway.) In Fig. [5] 
the error bars never do quite overlap at x/m ~ 0, but 
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they do overlap at x/m « ±1, although we must keep in 
mind that they are rough estimates. As discussed above, 
the error bars cannot be made equal at the origin by 
changing m/b, although the overlap can be made better 
by increasing the separation. The fact that the rescaled 
metric components in Fig. [6] take off in different direc¬ 
tions as they approach BH1 from between the holes is 
partly because the two metrics blow up at different co¬ 
ordinate locations. This small relative translation is due 
to the coordinate system used in the matching. 


B. Transition Functions 


The fact that matched curves do not strictly overlap 
even in the buffer zone means that a piecewise metric 
such as Eq. © possesses discontinuities at the transi¬ 
tion radii r A , wherever they are chosen to be. These dis¬ 
continuities can be problematic for numerical evolutions 
of the spacetime and thus it is desirable to smooth them. 
We now construct transition functions that smooth these 
discontinuities out, by letting 

9fc 3) = [l~F 1 (i?r)] g$ + F, (R 1 )g $, (41) 

where F\ has the properties that Fi(Ri > b) = 1, 
Fi(i?i < TOi) = 0, Fi(Ri « rf) ~ 1/2. This ansatz 
yields a metric that is equal to the inner zone metric 
near black hole 1, while it is equal to the near zone met¬ 
ric far away from black hole 1. In between (i.e. in the 
buffer zone) we obtain a weighted average of these two 
solutions. Since the Einstein equations are nonlinear, 
the sum of two solutions is in general not another solu¬ 
tion. But since both solutions are valid in the buffer zone, 
and since both have been matched, these two solutions 
are equal to each other in the buffer zone up to uncon¬ 
trolled remainders of 0(2, 3), corresponding to higher or¬ 
der post-Newtonian and tidal perturbation terms. Hence 
any weighted average of these two solutions in the buffer 
zone will yield the same correct solution up to uncon¬ 
trolled remainders of 0(2,3). This justifies the use of 
smoothing to merge the two solutions in the buffer zone. 
A similar solution can be obtained for the other black 
hole by replacing 1 —> 2. 

We choose transition functions of the form 


/ 0 ) = 


0 , r <r 0 

\ {i + tanh ii ( tan (^j0 - n>)) 


tan (^( r ~ r o)) , 

1 , r > rg + w. 


t)]} 


r 0 < r < ro + w 


(42) 


This function transitions from zero to one in a transition 
window which starts at rg and has a width w. The pa¬ 
rameter q controls the location at which / reaches 1/2, 
and s controls the slope at that location. Note that this 
transition function is C°°, a property which is useful for 
numerics. In Eq. m we set 


(43) 


g with transition function along the x-axis 

(tri | /m=0.5, b=10m) 



FIG. 7: In this figure we show the xx component of the global 
metric compared to the near-zone PN metric and the inner- 
zone BHPT metric around BH1. Observe that the transition 
function takes the global metric smoothly from one to the 
other. Error bars in the global metric are the same as the 
error bars of whichever local approximation is better at that 
point. 


with parameters 

ro = 1.5m, w = 5r A , g = 1/4, s = 10, (44) 

where r T A is given by Eq. m- With these parameters 
the transition functions reach the value 1/2 at Ra ~ 
1.5to + 0.16u> ~ 5.2m, very close to r A « 4.8m. We 
could have used something more sophisticated such as 
a transition with the same anisotropic behavior as the 
buffer zones, but several trials show that such details 
do not matter as long as the transition function has the 
right asymptotic properties. (Trials also showed that the 
results were not too sensitive to the precise parameter 
values.) 

We can see the effect of this transition function in Fig. [3 
and[8l In these figures we only show the region near BH 1, 
but similar behavior is observed near the other hole. The 
transition function effectively takes one solution into the 
other smoothly in the buffer zone. The size of the transi¬ 
tion window can be modified by changing the “thickness” 
w, but if w is made too small the derivatives of the metric 
(which include 1/w terms) develop artificial peaks inside 
the transition window. 

If the separation of the black holes is large enough 
so that the two transition windows of width w do not 
intersect, the transition functions given in Eq. (14211 will 
suffice to generate a global metric of the form 

9$ obal) = F 2 (R 2 )F 1 (R 1 )g$ + (l-F 1 (R 1 )\g$ 

+ [1 -F 2 (Ra)]g$, (45) 

However, if the separation is small enough that the two 
transition windows of width w start to overlap, we must 
construct a third function G(x) to allow for a smooth 


F a (R a ) = f(R A ) 
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g with transition function along the x-axis 

UU (m/n^O.5, b=l()m) 



(The more complicated form of ro is to account for the 
origin of the x coordinate being at the binary center of 
mass rather than on either hole.) 

The reader might worry that the use of transition func¬ 
tions could introduce large artificial gradients. However, 
with a reasonable choice of transition functions this is not 
the case. In order to understand why, let us look at the 
derivatives of the smoothed metric in more detail. Con¬ 
sider for example the overlap region O13 in which F 2 = 1 
and the smoothed metric is thus Figffl + [1 — Fi] g $. A 
derivative of this smoothed metric takes the form 


F i9 ^'+ [1-F 1 )g^' 


K 


o(3) 


< 7 (1) 

flu 


(49) 


FIG. 8: Same as the previous Figure, but for the 00 metric 
component. 


transition between the two black holes while not contam¬ 
inating the global solution near BH1 with a piece of the 
solution from inner zone 2 and vice versa. (In essence, 
this is an after the fact way of handling the fact that 
the buffer zones are not really spherically symmetric as 
often implied in the literature. When the holes are too 
close, a good transition function should not be spheri¬ 
cally symmetric either.) For the system considered here, 
b = 10 to is a sufficiently small separation that such a 
third transition function is necessary. The global metric 
then becomes 

9$° bal) = G(x) {F^R^J + [1 - FiORO] g$ } 

+ [1 - G(z)] {F 2 (R 2 )g$ + [1 

-F 2 (R 2 )]g$}. (46) 

(Recall that x is the distance along the axis between the 
holes with the origin at the center of mass.) The function 
G(x) will be chosen such that it is equal to unity near 
black hole 1 and zero near black hole 2. In between the 
two black holes, G{x) will range from zero to one, so 
that non-trivial averaging will occur only in this region. 
Again, this averaging is allowed because the two solutions 
in the curly brackets are both valid (and equal up to 
uncontrolled remainders) in between the two black holes. 
Near each black hole we obtain the appropriate inner 
zone solution, while far away Fi(I?i) = F 2 (I? 2 ) = 1, so 

(3) 

that we obtain the near zone solution gjjJ. We choose 
the transition function G(x) between C\ U C3 and C 2 U C3 
to be of the same form as the function in Eq. (14U1) . i.e. 

G(x) = /{x), (47) 


but with different parameter values 


r 0 = 


6 (to 2 — mi) b — m 
2 to 2 


, w = b — to, q = 1, s = 5/2. 

(48) 


The last term is the worrying one since it involves a 
derivative of the transition function. But note that 
the coefficient multiplying this term is the difference be¬ 
tween the inner-zone and near-zone metrics in the buffer 
zone, and therefore is by definition of the order of the 
uncontrolled remainders in the expansions of the first 
two terms. Therefore the third (unphysical) term will 
be safely absorbed into the small uncontrolled remain¬ 
ders unless we make a pathological choice of transition 
function—for example, one that has an inverse power of 
a small expansion parameter built into it. Our transition 
functions are chosen to avoid this. Their maximum slope 
is roughly l/rj, comparable to the slopes of g$ and gj£) 
in the vicinity of rj where the maximum occurs. Thus 
the unphysical third term in Eq. (1491) is always formally 
small. (We demonstrate that it is also small in practice 
in Sec. eh) The errors introduced into the constraint 
equations are found by differentiating ( 1451 ) . Again, all 
terms involving derivatives of the transition functions are 
multiplied by the differences between expanded metrics 
in the buffer zone, which are of the same order as the 
uncontrolled remainders and therefore can be neglected. 

In Figs. [9] and [TUI we demonstrate the good behavior 
of the smoothed solution. Observe that as x increases, 
the global metric with the transition function becomes 
identical to the near zone metric, while as ta approaches 
zero (near each hole, x/m —> ±5) it becomes equal to 
the inner zone metric. Since the holes are close to each 
other, in the region between the holes the global met¬ 
ric never becomes identical to the near zone metric but 
rather always contains a contribution from the inner zone 
metrics. This linear combination is valid in that region 
because there the asymptotic expansions of both approx¬ 
imate metrics are comparable to each other. In Figs. [7] 
and[H the error bars of the global metric with the transi¬ 
tion function overlap the error bars of the inner and near 
zone metrics in the regions where the former are valid. 
This criterion would not be satisfied if the near zone and 
inner zone curves were farther from each other, which 
would occur if the holes were closer, and thus could be 
taken as another indicator that the holes are still (barely) 
far enough apart for matching. 
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Global g along the x-axis 

(m^m=0.5, b=10m) 



FIG. 9: This figure shows the xx component of the global 
approximation of the metric across both holes. 


Global g along the x-axis 

(mj/m=0.5, b=l()m) 



FIG. 10: Same as the previous Figure, but for the 00 metric 
component. 


VII. INITIAL DATA FOR NUMERICAL 
RELATIVITY 

The approximate metric (1461) could be used as initial 
data for binary black hole simulations. To facilitate this 
task we now present this metric in the 3 + 1 decompo¬ 
sition, by providing explicit analytic expressions for the 
extrinsic curvature, lapse and shift on the T = t = 0 
slice. If the normal vector to this slice is denoted by n a , 
then the intrinsic metric in the slice is given by 


where £ n is the Lie derivative in the direction normal to 
the t = 0-slice. Below we compute K ab using the explicit 
expression 


£^ab — ^ 71 (&cQab &aQcb dbQac) 3 (52) 

which has been obtained using the ordinary derivative 
operator and the fact that q a bn b = 0. The evolution 
vector 


(d t ) a = an a + /3 a (53) 


is split into pieces perpendicular and parallel to the t — 0 
slice, where a denotes the lapse and (3 a the shift. Note 
that n a n a = —1 and /3 a n a = 0. 

The near zone extrinsic curvature computed from the 
PN metric [3lj in corotating harmonic coordinates on the 
t = 0 slice is given by 


A+ = £ 

A =1 


*7 


m A 


^ V A n A - ^ 


(54) 


where the error is of 0(m/b) 5 / 2 , the superscript 3 is to re¬ 
mind us that this expression is only valid in C3, and where 
the parentheses on the indices stand for symmetrization. 
In the previous equation, v A and n l A denotes the particle 
velocities and directional vectors given by 


2 m 2 

Ui = ui -0, 


2 

v 2 =-u — b, 
m 


v\=v\ = 0 (55) 


and 


k _ x X A 
n A — I 

ta 


1 TO 2, 1 

X 1 = - 0 , x 2 = 


— — -b, x A =x A =0. 
m 


The corresponding near zone lapse and shift on the t = 0 
slice are 


«(3) 


2 


i-E 

A=1 


mA 

XA 


(57) 


and 


ai ^m A v l A a 

P(3) = - 2 ^ ~r. - ewux 3 


A =1 


ta 


(58) 


qab = gab + n a n b , (50) 

and the extrinsic curvature is 

-fc-ab — ~z£nqabi 


where once more this is valid on C3. 

The extrinsic curvature of the T = 0 slice valid in the 
inner zone of black hole 1 (Ci) and computed from the 
metric given in the previous section in isotropic corotat¬ 
ing coordinates is 







19 


K, 


(1 ),ICC 
00 


m2 


= --ptfYiX'+Y*) 


4/5 


K, 


(1 ),ICC 
01 


* M 3nx 

b 4/ _ Ml 
2 + Y 2 - Z 2 ) 


4/ _ Mi. 

* Ri 




Y 2 f, A 2Ml 


b 2 


- 


R 2 


+ Hk- 


Mi 

~rI 


^X 2 + 2Y 2 )+ U-^) Ml ^ b 


K, 


(1 ),ICC 
02 


(X 

m 2 XY 

b 3 


n 


4> 5 


4/ _ Ml 

* fli 


(* - £)'^ + 3 “ (* 4 - 2 ~w) - (* - fr) w“ <** + ?2 - 22) 


A^ )JCC = ^nxzv* 


Mi /—9 —9 —9 \ M\ 

— (x 2 + y 2 ~z 2 )-* + — 


r/ (i),;cc _ m' 2 ^ 'J' 5 

A H “ OJ.3 vjy 

U'hX 2 


2 b 3 H> - Ml 

ft-i 


4 ME' - 


Mi 

in 


bu + 4 4> - 


M A.Mi ^ o,t,4 , Mix 2 




+ - 


R\ 


K 


(1 ),icc _ m 2 X 41 5 


12 


26 3 

ft 1 

(x 2 - y 2 - z 2 


M \ 2 XY 2 

2(4- -- ) 6w + 12Mift—=- 

y i?i 


i?? 


1 + 


wl" - 3HA 4 


i? 2 


- 24' 4 + 




M 2 


. T Mi \ , Mi 

2[V& ^J 6 ^ 


r (D,/cc 6Mim 2 flX 2 YZ T 5 / M 2 \ 

13 6 3 i?f ml \ + 4R 2 J ’ 


K. 


(1),/CC 


m 2 T ik 5 


K. 


(1),ICC 

23 


v 33 


26 3 4- 

Mi ) 

Ri ( 

m 2 Z 

V]/ 5 

2b 3 


m 2 Y 

4/5 


12^-HA - 6J]f 4 X + 3 Ml3 f ELY 

2(5-- ^4 bu - 12Mifi^^ 

1 M ) i?? 


i?? 


3HAA 2 


$ - — ) buj (A 2 - Z 2 ) 

i?i ' 


(59) 


1 + ^- 2 ^-^^(^ + F2 ^ 2 ) 


i?l / Rf 


263 

Ml 


4 ME' - 


Mi 

in 


Mi 


buJ+ 41 * - buj ^~ + ( 4 if - 24,4 


MfZ 2 M\Z 2 

+ 4 




R\ 


r 


where the error is of 0(R±/b) 3 , and where the superscript 
ICC is to remind that that this is calculated in isotropic 
corotating coordinates. Later on, we will transform this 
metric to harmonic corotating coordinates with the map 
(f>i 3 found in Eq. m , and we will drop this superscript. 
In the above equations, 4/ is the Brill-Lindquist factor for 
black hole 1 in isotropic coordinates, i.e. 


coordinate transformation and are obtained from the 
purely spatial components K^’ ICC using 


t A1),igc _ k „l K {1 ),ICC 
/v 0i/ — % c lv 1 '-kl > 


where the projection tensor q” is given by 


(61) 


T = 1 


Mi 

2 Rl' 


(60) 


The K^’ ICC components will be needed later in the 


< = ^ + n ^ )JCCn (i) ; icc- ( 62 ) 

Here, the normal vector to the slice computed with 
(1 ),ICC 


= -\/- 1 /h° ( l ) (dT) a is given by 


a),icc 


1 _ Mi. 

2Ri 

T 




nP’ ICC = 0, 


(63) 


where P 2 stands for the second Legendre polynomial. 
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The upper components are then 
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(64) 


This means that the lapse and shift of the T = 0 slice in 
the inner zone (C i) and in inner corotating coordinates 
are given by 


Q (1),/CC — ~ n 0 


(1 ),ICC 


I _ Ml 

2fli 

T 






/cc 


0^p 2( x /R i) 

- a (i),iccri(i)jcc- 


(65) 


Observe that the lapse of Eq. (l65l) goes through zero at 
Pi = M/2. Apart from a small perturbation it closely 
resembles the standard Schwarzschild lapse in isotropic 
coordinates. 

Note that the extrinsic curvature, lapse and shift given 
up to this point are expressed in two different coordinate 
systems. The post-Newtonian quantities valid in the near 
zone (C 3 ) are given in harmonic corotating coordinates, 
while the black hole perturbation theory results valid in 
the inner zones (Ci >2 ) are given in isotropic corotating 
coordinates. We will now apply the coordinate transfor¬ 
mation found in Sec. El namely Eq. (ED, to transform 
the inner zone expressions to harmonic corotating coor¬ 
dinates, thus dropping the label ICC in favor of the su¬ 
perscript (1). The result for the inner extrinsic curvature 
of black hole 1 is given by 




( 66 ) 


where all components still have errors of 0(5/2,3). The 
extrinsic curvature in harmonic corotating coordinates in 
submanifold C 2 can be obtained from the above equation 
by the symmetry transformation discussed in Eq. (1321) . In 
Fig. HD we have plotted the xy component of the extrinsic 
curvature along the a;-axis. Observe that the behavior 
of the post-Newtonian solution (dashed line) is different 
from that of the black hole perturbation solution close to 
the black hole, where the latter diverges more abruptly. 
In Figs. EHISl the error bars have been estimated by 
inserting Eqs. ED and (138|) into the definition of the 
quantity plotted and using dtrA ~ bio and dtRA ~ bfl. 



x/m 


FIG. 11: This figure shows the xy component of the near zone 
(PN - dashed line) extrinsic curvature, as well as the inner 
zone curvatures (BHPT - dotted lines) obtained via black hole 
perturbation theory. This figure uses the same test case as 
previous figures, with equal mass black holes and m/b = 1/10. 


Similarly, the lapse and shift in harmonic corotating 
coordinates corresponding to the inner zone of black hole 
1 are given by 

1 (1) t T 

“(1) - Zo~ - ~ n 0 - J t(X{l),ICC , 

n (l) 

P[ 1 ) = -ami, (67) 

where again the lapse and shift for the inner zone around 
black hole 2 (C 2 ) can be obtained by the symmetry trans¬ 
formation (1521) . In these equations, the normal vector is 
given by 


n (i) — J” v n \i )jcci (6®) 

where the matrix J^ v has been defined in Eq. (1551) . Note 
that «(!) in Eq. (1671) has the same zeros as «(i ),icc an d, 
thus, «(i) also changes sign at R\ = M/2. Furthermore, 
since J T t = 1 + 0{m/b) the inner zone lapse a(i) equals 
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x/m 


FIG. 12: This figure shows the near zone (PN - dashed line) 
and the inner zone lapse (BHPT - dotted line) along the x- 
axis. Observe that the near zone lapse crosses zero at x/m = 
5.5 and x/m = 4.5, which is the location of the event horizon 
in harmonic coordinates. 


a (i),icc U P to a perturbation of 0(m/b) and thus a(i) is 
equal to the standard lapse of Schwarzschild in isotropic 
coordinates plus a perturbation of 0(m/b). These fea¬ 
tures are borne out by the plot in Fig. I~HZ1 which shows 
the global lapse along the x-axis. We can also see from 
the figure that while the inner zone lapse goes to —1 as 
7T —> 0, the near zone one diverges to negative infinity. 

With these equations, we can construct an approxi¬ 
mate piecewise global extrinsic curvature, lapse and shift 
by substituting the metric for these quantities in Eq. (1361) . 
By the theorems of asymptotic matching, the derivatives 
of adjacent pieces of the piecewise metric will be asymp¬ 
totic to each other inside their respective buffer zones. 
This asymptotic similarity is, thus, also observed in the 
extrinsic curvature, as well as the lapse and the shift. 
Due to the piecewise nature of these solutions, there will 
be discontinuities on a 2-sphere located at some transi¬ 
tion radius inside of the buffer zone. In order to elim¬ 
inate these discontinuities, we use the same transition 
functions used for the metric in Eq. (|42|f with the same 
parameters. In this manner, we obtain a smooth global 
extrinsic curvature, lapse and shift given by 


J 


global) = G{x) { Fl{Rl)K if + [! _ Fl ( Rl)] K \f) + [1 - G(x)} { F 2 (R 2 )K + [1 - F 2 (R 2 )\ K[f) , 
oi(global) = G(x) {Fi(i?i)o!( 3 ) + [1 - -Fi(-Ri)] a(i)} + [1 - G(x)\ {F 2 (R 2 )a (3) + [1 - F 2 (R 2 )\ 0 ( 2 )} , 

^global) = G{x) [f^Ri)^ + [1 - F 1 {R 1 )\ /Jfo} + [1 - G{x)\ {f 2 {R 2 )(5\ z) + [1 - F 2 {R 2 )\ $ 2) } . (69) 


Figs. ii and M show the global lapse, shift and ex¬ 
trinsic curvature with the transition functions. 

We also could have computed the extrinsic curvature 
directly from Eq. This would add terms involving 

derivatives of the transition functions. The parameters of 
the transition functions were chosen so that these deriva¬ 
tives are of the same order as the uncontrolled remainders 
in the buffer zone, and thus formally do not affect the ac¬ 
curacy of the extrinsic curvature. Since the two methods 
are equivalent, we took the one which was simpler to 
compute (had fewer terms). 

The initial data constructed by the methods above 
[Eqs. (HSll and (1551) ] is only an approximate solution to the 
Einstein equations. Therefore, this data leads to an error 
in the constraints of the full theory of order 0(ni/b ) 9//2 
near the horizons and 0(m/r) 2 in the near zone. This 
error can be sufficiently small compared to other sources 
of numerical error such that solving the constraints more 
accurately is not required. However, perhaps the optimal 
approach would be to use this solution as input to York’s 
conformal method jhj and compute a numerical solution 
to the full constraints. Since this data is already signif¬ 
icantly close to the constraint hypersurface, there might 


be some hope that appropriate projection methods will 
not alter much the astrophysical content of the initial 
data. Somewhat surprisingly, standard PN data (with¬ 
out matching) has not yet been used for the generation 
of numerical black hole initial data except in Ref. El: 
which is based on the PN data of Ref. [33|. We leave it 
to future work to explore similar techniques for the data 
set presented here. 


If this data is to be evolved, it is necessary to choose a 
lapse and a shift. The choice presented in this section is 
natural in the sense that it is close to quasi-equilibrium. 
In other words, with the lapse and shift presented in 
this section, the 3-metric and extrinsic curvature should 
evolve slowly. However, since our lapse is not everywhere 
positive, some evolution codes may have trouble evolv¬ 
ing with it. If this is the case, one can simply replace 
the above lapse with a positive function at the cost of 
losing manifest quasi-equilibrium, but without changing 
the physical content of the initial data or the results of 
the evolution. 
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a along the x-axis 

(mj/mSO.5, b=10m) 



FIG. 13: This figure shows the global lapse along the x-axis 
with the transition function. 



FIG. 14: This figure shows the global shift vector along the 
x-axis with the transition function. 


VIII. CONCLUSIONS 

We have constructed initial data for binary black hole 
evolutions by calculating a uniform global approxima¬ 
tion to the spacetime via asymptotic matching of lo¬ 
cally good approximations. The manifold was first di¬ 
vided into three submanifolds: two inner zones (one for 
each hole), C\ and C 2 (+4 b) equipped with isotropic 
coordinates; and one near zone, C 3 (ta tua and 
r <C A/27r), equipped with harmonic coordinates. In 
the near zone, the metric was approximated with a post- 
Newtonian expansion, while in each inner zone the met¬ 
ric was approximated with a perturbative tidal expan¬ 
sion of Schwarzschild geometry. Each approximate solu¬ 
tion depends on small parameters locally defined on each 
submanifold. These submanifolds overlap in two buffer 
zones, 0 i 3 and 0 2 3 (4—volumes), given by the intersec¬ 
tion of the inner zones with the near zone, i.e. on an 


Global K along the x-axis 

(m^m^.5, b=10m) 



FIG. 15: This figure shows the global xy component of the 
extrinsic curvature along the x-axis with the transition func¬ 
tion. Note that the bumps due to the transition function 
are comparable to the error bars which estimate uncontrolled 
remainders in the expansions. 


initial spatial hypersurface the buffer zones becomes 3- 
volumes given by to a <C ta b. Inside each buffer zone, 
two different approximations for the metric were simul¬ 
taneously valid and hence we were allowed to asymptot¬ 
ically match them inside this zone. 

The matching procedure consisted of first expanding 
both adjacent approximate metrics asymptotically inside 
the buffer zones. After transforming to the same gauge, 
these asymptotic expansions were then set asymptotic to 
each other—equating their expansion coefficients, which 
does not in general set the functions equal to each other 
anywhere in the buffer zone. After solving the differential 
systems given by equating expansion coefficients, asymp¬ 
totic matching returned a coordinate transformation (^13 
and 4 > 23 ) between submanifolds and matching conditions 
(^13 and '0 2 3 ) that relate parameters native to different 
charts. A piecewise global metric was then obtained by 
transforming all metrics with the set {<t> n m', V’nm}, result¬ 
ing in coordinates which resemble harmonic coordinates 
in the near zone but isotropic coordinates in each inner 
zone. 

Once a piecewise global metric was found, the spatial 
metric and extrinsic curvature were calculated in each 
zone by choosing a spatial hypersurface, with the stan¬ 
dard 3 + 1 decomposition. This initial data was then 
transformed in the same manner as the 4 metric. Due 
to the inherent piecewise nature of asymptotic match¬ 
ing, this data was found to have discontinuities of or¬ 
der 0(3/2,3) or smaller inside the buffer zone. We con¬ 
structed transition functions to remove the remaining dis¬ 
continuities in metric components and spikes in deriva¬ 
tives. These transition functions were carefully built to 
avoid introducing errors larger than the uncontrolled re¬ 
mainders of the approximations in the buffer zones. With 
these functions, we constructed a global uniform approx- 
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imation to the metric valid everywhere in the manifold 
with errors 0(m/b ) 9 / 2 near the black holes and 0(m/r ) 2 
far away from either of them. 

This uniform global approximation of the metric can be 
used as long as the black holes are sufficiently far apart. 
When the two black holes are too close, there is no in¬ 
tervening post-Newtonian near zone in which to match. 
However, since there is no precise knowledge of the region 
of convergence of the PN series, it is unknown precisely 
at what separation the near zone vanishes. We have ex¬ 
perimented with separations b > 10m and we have found 
that, in these cases, a region does exist between the holes 
where the post-Newtonian metric is reasonably close to 
the perturbed black hole metrics and thus matching is 
still possible. For separations of b < 10m, this region 
shrinks rapidly and matching is not guaranteed to be 
successful. Also, the “global” metric is not valid all the 
way to the asymptotically flat ends inside the holes, im¬ 
plying that our initial data must be evolved with excision 
techniques rather than punctures. 

We then constructed a lapse, shift, and extrinsic cur¬ 
vature, all of which are needed for numerical evolutions. 
The lapse was found to possess the expected feature that 
it becomes negative inside the horizon of either black 
hole. Some numerical codes might find this feature un¬ 
desirable, in which case the lapse can be replaced by some 
positive function at the cost of losing approximate quasi- 
stationarity. These 3 + 1 quantities were then smoothed 
with transition functions of the same type as those used 
in the 4-metric. 

In conclusion, we have constructed initial data for an 
inspiraling black hole binary that satisfies the constraints 
to order 0[{m/b) 5 ^ 2 (Ri/b) 2 ] in the inner zone, and to or¬ 


der 0(m/r) 2 in the near zone. This data is a concrete 
step toward using PN and perturbation methods to con¬ 
struct such initial data, and it should be compared to 
other numerical methods with respect to its ability to 
approximate the astrophysical situation. We should note 
that the data presented here makes use of perturbative 
expansions of low order ( e.g ., the near zone metric is 
built from a 1 PN expansion), but this paper firms up a 
method introduced by Alvi [2C| that could be repeated 
to higher order at the cost of more algebra |T(| [13 • The 
post-Newtonian metric needed for the next order in m/r 
(and beyond) is available f3jj, as is the octopole pertur¬ 
bation ( r/b ) 3 of a Schwarzschild black hole (24[. Our 
method might also be extended to spinning black holes, 
which are more astrophysically realistic^ with the avail¬ 
able post-Newtonian near-zone metric 1301] and tidal per¬ 
turbation dj. 
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